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Abstract We outline the basic ideas and techniques underpinning the simulation of 
stochastic differential equations. In particular we focus on strong simulation and its 
context. We also provide illustratory examples and sample matlab algorithms for the 
reader to use and follow. Our target audience is advanced undergraduate and graduate 
students interested in learning about simulating stochastic differential equations. We 
try to address the FAQs we have encountered. 
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1 Introduction 

1.1 When is a model stochastic? 

Often in modelling, we need to incorporate a phenomenon or influence that seems ran- 
dom, whose behaviour we can only recognize as statistically Gaussian. The prototypical 
example is behaviour of molecular origin — the Brownian motion of a pollen particle 
on a water surface. However the influences can be large scale, for example a turbulent 
wind or atmospheric flow, or thousands of people buying and selling millions of shares. 

Imagine tracking and adjusting the flight of a rocket after lift-off. If the rocket 
is buffeted by a random turbulent wind, you might sensibly equip the rocket with 
stabilizers that kick-in if a gust diverts it too far. Computing the path of the rocket, and 
regulating it to ensure it threads pre-arranged target positions at critical junctures (eg. 
stage separation), is a stochastic simulation problem. Indeed it is a strong simulation 
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problem, as conditional on the path, the stabihzers will affect the outcome. The pricing 
of financial derivatives (futures/options) is another example. One tracks a security price 
(eg. a share price) that is randomly buffeted by market forces. Pricing the derivative 
you wish to sell, which might be exercised by the buyer at a future time, involves 
hedging/regulating the proportion of your investment in the security (the rest invested 
in risk-free bonds) so as to minimize your risk. Building in stabilizers/barriers that 
kick-in if the security price skews too wildly, is again, a strong stochastic simulation 
problem. 



1.2 What is a stochastic difTerentiai equation? 

Consider a model that incorporates some random phenomena whose statistics is Gaus- 
sian. Suppose the state of the system is recorded through the vector yt G K^, with 
N 1. Suppose there are several random sources, say W'^,...,W'^; these are Wiener 
processes; think of them as independent, continuous, nowhere differentiable, functions 
of time. Indeed the time derivatives of the Wiener processes represent pure white noise. 
Suppose the affect of the Wiener processes on the model, i.e. which way they skew the 
solution, is recorded through the vector fields Vi, . . . , V^. Without the noise we would 
have a nice fuzz-free signal which is generated by the vector field Vq. A stochastic 
model for the system state yt € might be that it evolves according to the stochas- 
tic differential equation: 

This representation of the model is somewhat formal; after all the pure white noise 
terms dWl /dt need to interpreted in an extremely weak sense, we prefer to represent 
the model in the form 

dyt = Voiyt) dt + Vi{yt) dWl + ■ ■ ■ + Vaiyt) dwf. 

Indeed there is also a representation in a preferred integral form which we meet 
presently. With this in mind though, an important observation at this point, is to 
recall that Wiener processes are continuous functions. Thus the solution function will 
also be continuous. 

Example (Langevin equation/Brownian motion). Consider the equation of motion of a 
pollen particle suspended in a fluid flow. The particle might obey the following equation 
of motion for its velocity yt: 

dyt , rr dWt 

_^-ayt + Vb—, 



where a and b are constants. The right-hand side is the force exerted on the particle 
per unit mass. There is a deterministic force —a yt and a white noise force VbdWt/dt 
supposed to represent the random buffeting by the water molecules. 
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1.3 What information do we want to retrieve? 

If the time interval of interest is [0,T], and our initial deterministic state is yo, each 
realization uj of an individual Wiener path W{lj) will produce a different outcome yi(ui) 
for t £ [0, T]. Practical information of interest is often expected values of functions 
/ of the solution, f{yt), or more generally path-dependent functions of the solution 
f{t, ys;s ^t). Hence we might want to compute 

E/(2/t) := J f{yt{oj)) dPH, 

where P is a probability measure. For example we could pick / to be of polynomial or 
exponential form, synomynous with statistical moments of yt- If / is the identity map, 
we obtain the expectation of the solution. If we take / to be || • \\p, where || • \\p is the 
p-vector norm, then we define the L^-norm for p ^ 1 by 

\\yt\\% ■■= j 112/tHii^dPH. 

1.4 hiow do we retrieve it? 

There are two main simulation approaches to extract such information, we can either: 

— Solve a partial differential equation-^ or 

— Perform Monte-Carlo simulation. 

Associated with every stochastic differential equation, there is a parabolic partial 
differential equation for u{t, y) whose solution at time t € [0, T] is 

u{t,y) = Ef{yt) 

provided u(0, y) = f{y) initially. Thus solving the associated partial differential equa- 
tion on [0, T] will generate lots of information about the solution to the stochastic 
differential equation at time t. By dudiciously choosing / to be a monomial function 
we can generate any individual moment of the solution yt we like, or if we choose 
/ = exp we generate all the moments simultaneously (this is essentially the Laplace 
transform). If we choose / to be a Dirac delta function we generate the transition 
probability distribution for yt — the probability density function for yt conditioned on 
the initial data j/o- Choosing / to be a Heaviside function generates the corresponding 
(cumulative) distribution function. Of course often, the partial differential equation will 
have to be solved approximately. Also note, if we fix a form for / from the start, for 
example / as the identity map, then we simply solve an ordinary differential equation 
for u{t,y). 

In Monte-Carlo simulation, we generate a set of suitable multidimensional sample 

paths W(iLi) ■— {W^{u)), . . . ,W'^{uj)^ on [0, T]; in practice, uo belongs to a large but 
finite set. For each sample path W{uj), we generate a sample path solution y{uj) to the 
stochastic differential equation on [0, T]. This is often achieved using a truncation of 
the 'stochastic' Taylor series expansion for the solution y of the stochastic differential 
equation, on successive small subintervals of [0,T]. Suppose for example, we wanted 
to compute the expectation Efijjt). Having generated a set of approximate solutions 
yti'jJi) at time t £ [0, T], for i = 1, ...,P with P large, we can estimate Ef{yt) by 
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computing the mean-sum over the large finite set of approximate sample solutions 
yt{y>i). Hence in practice we approximate 



/ /(j/t(a;))dPH^^^/(ytK)) 

where P is the total number of sample paths. A natural dichotomy now arises. To 
compute Ef{yt), we can in fact choose any suitable multidimensional paths W{lj) that 
leave E/(yt) approximately invariant, in the sense that \\E f(yt) — E f{yt)\\ is sufficiently 
small. This is a weak approximation. For example, increments AW^ in each computa- 
tion interval can be chosen from a suitable binomial branching process, or using Lyons 
and Victoir's cubature method [28]. Note that since the approximate paths are not 
close to Brownian paths we cannot compare yt and yt directly. In a strong approxi- 
mation, discrete increments AW^ in each computation interval are directly sampled 
from the Gaussian distribution. This is more expensive. However, the sample paths 
W{uj) generated in this way, allow us to compare yt(a;) and yt(i^) directly in the sense 
that we can guarantee E ||yt — yt\\ will be sufficiently small. Naturally, using strong 
simulation we can also account for path-dependent features, such as conditional cut- 
offs or barriers, when we investigate individual solutions or the final expectation or 
higher moments of the approximate solution paths y. For a comprehensive overview of 
Monte-Carlo methods see Boyle, Broadie and Glasserman [6]. 



1.5 What is required? 

In general to extract qualitative and quantative information from a stochastic dif- 
ferental system requires the languages and techniques of several mathematical disci- 
plines, notably: 

1. Integration: in Brownian motion new information is continuously generated on in- 
finitesimally small time scales (imagine the pollen particle jiggles); solution as with 
ordinary differential equations is by integration, except that now the coefficients of 
the evolution equation — the Wiener processes — are no longer differentiable. 

2. Statistics: we typically extract statistical information from the solution process; 

3. Geometry: as with ordinary differential equations, preserving invariant geometric 
structure of the solution path evolution is important; for example the solution may 
evolve on a homogeneous manifold; 

4. Simulation: stochastic differential equations more often than not, are not integrable 
in the classical sense, and require numerical computation. 

For general background reading, we recommend as follows. For a comprehensive in- 
troduction to the theory underlying stochastic differential equations download Evans' 
notes JTSJ. For an introduction to numerical simulation, see Higham's notes [19]. The 
answer to just about any other question that a beginner may have on numerical sim- 
ulation, not covered above or here, can likely be found in the treatise by Kloeden and 
Platen 122. 
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2 Stochastic differential equations 

2.1 Integral representation 

Consider the nonlinear stochastic differential equation of order N £ N given by 

/■t d .t 

yt = yo+ / Vb(yr)dr + V / v^iyr)dWl 
Jo Jo 

Here {W^,...,W'^) is a d-dimensional Wiener process, i.e. there are d independent 
driving noisy signals. We assume there exists a unique solution J/: [0,T] i-^ for some 
time interval [0, T] C R+. We suppose that Vq and : , i = 1, . . . , d, are 

smooth non-commuting autonomous vector fields. We are representing the stochastic 
differential equation above in Ito form and indicate this by using Vq to represent the 
ltd drift vector field. We call the vector fields Vi for i — 1, . . . ,d associated with the 
driving noise terms the diffusion vector fields. Presently we will distinguish, and explain, 
the Ito representation as opposed to the Stratonovich representation for a stochastic 
differential equation. We also remark that a common convention is to set = t. 
Results on existence and uniquess of solutions can be found in Kloeden and Platen [22] . 

2.2 Driving Wiener process 

A scalar driving noisy signal or disturbing Brownian motion has a concise definition 
and set of properties formulated by Wiener. 

Definition 1 (Wiener process) A scalar standard Wiener process or standard Brownian 
motion Vl^ is a continuous process that satisfies the three conditions: 

1. Wq = with probability one; 

2. Wt-Ws ~ - s • N(0, 1) for < s < t, where N(0, 1) denotes a standard Normal 
random variable; 

3. Increments Wt — Ws and — Wrj on distinct time intervals are independent, i.e. 
ioT ^ s < t < rj < ^. 

Note that with probability one an individual Brownian path is nowhere differentiable. 

Example (Langevin equation). As we have seen, the Brownian motion of a pollen particle 
suspended in a fluid flow obeys the following equation of motion for its velocity yt : 

dyt = -ayt dt + VbdWt, 

where a and 6 are constants, and W is a scalar Wiener process. This type of stochas- 
tic differential equation is said to have additive noise as the diffusion vector field is 
constant. It is also an example of an Ornstein-Uhlenbeck process. 

Example (scalar linear equation). Consider the scalar linear stochastic differential equa- 
tion 

dyt = aytdt + byt dWt 

driven by a scalar Wiener process W, with a and b constants. This stochastic differen- 
tial equation is said to have multiplicative noise as the diffusion vector field depends 
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(N-1)h T=Nh 



Fig. 1 Example scalar Wiener path on the interval [0,T]. 



multiplicatively on the solution yt- We can in fact analytically solve this equation, the 
solution is 

yt = j/o exp(af + bWt- \b^t) . 
The additional term ' — ib^t' is due to the Ito correction, which we discuss shortly. 



2.3 Vector fields and flow maps 

Consider an ordinary differential equation governed by an autonomous vector field V 
that evolves on a homogeneous manifold M so that 



Ayt 
At 



= V{yt), 



with initial data yo £ M. Here we suppose M to be an A'^-dimensional manifold. The 
reader can for simplicity assume M = for the rest of this section if they choose. Let 
Diff(AI) denote the group of diffeomorphisms of A4. The flow-map ^pt,to £ Diff(A^) for 
the ordinary differential equation above is the map taking the solution configuration 
on M at time tp to that at time t, i.e. it is the map ipt,to '■ M ^ M such that 



'pt,to ■ yto 



yt- 



In other words for any data yo £ M ai time to we can determine its value yt £ Ml at 
time t later by applying the action of the flow map ipt,to to yo so that yt — ft.to ° J/O- 
Note that the flow map satisfles the usual group properties 

^t,s ° •^s,to = Vt,to> 

with 1(5*0, to = id, the identity diffeomorphism. If / G Diff(A^), the chain rule reveals 
that for all 2/ G we have 



d/(y) 
dt 



V{y)-dyf{y), 
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where dy = Vy is the usual gradient operator with respect to each component of y. In 
other words, vector fields act on the group of diffeomorphisms Diff(A4) as first order 
partial differential operators, and for any / e Diff(AI), we write 

Vofoy = V{y)-dyf{y). 

We now think of V{y) ■ dy as a first order partial differential operator and an element 
of the tangent space X(M) to T)iS{M). In particular, choose the diffeomorphic map / 
to be the flow map ipt = ftfl- Then for all yo € M and yt = ° yo we have, 

^ift oyo) = V oipto yo. 

We pull back this ordinary differential equation for yt € M. to the linear functional 
differential equation in Diff(AI): 

d</?t 

Since (po = id, the solution is (ft = exp{tV), giving the represention of the flow-map 
as the exponentiation of the vector field. Hence we see that 

yt = exp{tV) o yQ. 

An important and illustrative concomitant derivation of this result is as follows. 
Integrating the functional differential equation for the flow-map we get 

/■* 

(ft = id+ V o (fiT dr. 
Jo 

To solve this integral equation, we set up the formal iterative procedure given by 

(fil — id + / V o ip). ^ dr, 
Jo 

(0) (2) 12 2 

with (fl ' = id. For example, after two iterations: (p). ' = id + tV o id + ^t V o id, 

where V'^ = V oV . Hence in the limit we obtain the exponential form for '^i above. 
The composition of two vector fields U and \^ is a second order differential operator: 

U<^V={U{y)-dy){V{y)-dy) 

= ((f/(2/) • dy) {V{y))) ■ dy + U{y) ® V{y) : dyy 

N N 

= U'dy,{V^)dy^ + J2 U'V^dy.y,. 

Importantly we now observe that since dy^y^ = dy^y^ as operators on Diff(AI), the Lie 
bracket [U, V] of two vector fields is a vector field: 

[U,V] = U oV -V oU = (^{U{y) ■ dy)V{y) - (Viy) ■ dy)U{y)y dy. 

Note the sum of two vector fields is itself a vector field. Hence the set of vector fields 
is a Lie algebra X{M.) — closed under summation and Lie product [•,•]. 
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2.4 Stratonovich representation 

There are two generic representations for stochastic differential equations. One can 
either express them in Ito form, as we did at the beginning of Section 2, or we can 
express the stochastic differential equation in Stratonovich form, in which case we write 

Jo .^-^Jo 

Two subtle notational changes can be spotted. First the ' " dW^' indicates that the 
stochastic integrals on the right are supposed to be interpreted in the Stratonovich 
sense; discussed presently. Second we now use the Stratonovich drift vector field Vo 
instead of the Ito drift vector field Vq. The relation between the two vector fields is 

d 

Importantly, when stochastic integrals are interpreted in the Ito sense, then they are 
limit of a left Riemann sum and when repeated integrals are computed, an Ito cor- 
rection must be taken into account; a practical discussion of this point can be found 
in Higham |19l pp. 530-1]. For example, the correct evaluation of the an Ito repeated 
integral of a Wiener process with respect to itself is 

/ WldWl^\[wirf -\t. 
Jo 

Stratonovich integrals are interpreted as the limit of the midpoint rule, so for the 
corresponding Stratonovich integral the correct evaluation is 

/ W^.''dW^ = ^{wi^)^. 
Jo 

Thus the rules of Stratonovich integral calculus match those of standard integral cal- 
culus. For this reason it is often preferable to use the Stratonovich representation for a 
stochastic differential equation. The two representations are equivalent, but it is impor- 
tant to know in which form you have been quoted the stochastic differential equation. 
Often this depends on the modeller and their field of interest. In finance applications 
the Ito representation predominates, and by simply replacing the given Ito drift vector 
field Vq by the corresponding Stratonovich drift vector field Vq above, one can pro- 
ceed using standard integral calculus rules. In physical applications, the model is often 
often directly expressed in Stratonovich form. Hereafter we will use the Stratonovich 
representation and omit the ' ° ' symbol, unless we specify otherwise. 

3 Stochastic Taylor expansion 

We follow the procedure we performed above for the ordinary differential equation to 
try to find the solution for the fiow-map. In the process we obtain a solution series 
expansion called the stochastic Taylor series. We define the flow-map tpt £ Diff(AI) for 
the stochastic differential equation above as the map taking the solution configuration 
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on Ai at time to that at time t: hence yt = Vo- Using the Stratonovich represen- 
tation for a stochastic differential equation and the convention = t, the chain rule 
for any function / e Diff(AI) yields the stochastic differential equation governing the 
evolution of / o j/j as follows 

rt 

/o2/t = /oj/o + V / {Vi-dyf)oyrdWi 

As for the ordinary differential equation, setting f — ipt, we can pull back the stochastic 
flow on AI to a functional stochastic differential equation on Diff(AI) given by 

d ,.t 

V3f = id + V / ViOifir AWr- 
To solve this equation, we set up the formal iterative procedure given by 

with iff^ = id. By performing the iterations one can see formally, and prove rigorously, 
that the solution flow-map is given by the series expansion 

= id + y^{wi) Vi+y^( f r AWi^ AW^^ + • • • . 

i=0 i,j=o^^o ^0 ^ 

Here we use the notation Vij = VioVj. We can apply this to the initial data yo € M 
and obtain the stochastic Taylor expansion for the solution 



yt = yo + y^(V^t) Viiyo) +y,( I r aw;, AW^) Vij{yo) + ■ 



We can express the solution series for the flow-map concisely as follows. Let A* denote 
the free monoid of words over the alphabet A = {0, 1, ... , d}. We adopt the standard 
notation for Stratonovich integrals, ii w = ai . . . an then we set 



Jw{t) ■■= [ •••/"' AW?^ ■ ■ ■ AW?,". 
Jo Jo 



We also write the composition of the vector flelds as Vw = Voi o Va, o • • • o Va„ ■ Then 
the flow-map is given by 

ft = ^ Jw{t) Vw 
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4 PDE simulation 

There is an intimate link between any stochastic difTerential equation and a prescribed 
parabolic partial differential equation. The link is given by the Feynman-Kac formula, 
which we give here in a very simple form. See for example Karlin and Taylor [211 
pp. 222-4] for the full statement of the Feynman-Kac formula and its applications. 

Theorem 1 (Feynman-Kac formula) Consider the parabolic partial differential equa- 
tion fort G [0,r]; 

dtu = Cu, 

with u{Q, y) — f{y). Here £ := Vq + \ {Vi + • ■ ■ + V^) is a differential operator of order 
2N . Let yt denote the solution to the stochastic differential equation for t G [0, T]: 

yt=yo+ / Vo{yr)dT + y] / V^{yr)dWl 
Jo ,^1^0 

Then, when yQ — y we have: u{t,y) — Ef{yt). 



Remark. Note that using the relation between the Ito and Stratonovich drift vector 
fields, an equivalent formulation is C = Vq ■ dy + ^ ® '^i) ■ ^vv- 

We provide a purely combinatorial proof of the Feynman-Kac formula. Before we 
begin, we need the following results, for the expectation of Stratonovich integrals and 
also a combinatorial expansion. Let D* C A* denote the free monoid of words con- 
structed from the alphabet D = {0, 11, 22, . . . , dd}. The expectation of a Stratonovich 
integral Jw is given by 



EJu 



2<i(™)n(«))! ' 

0, 



w £ D"; 
w e A*\I 



In the formula, A{w) is the number of non-zero consecutive pairs from D in w and 
n{w) = z(w) -)-d(w), where z(w) is the number of zeros in w. We also have the following 
combinatorial identity for all G D*: 



{Vo + \{v^ + --- + vhf^ J2 



n{w)—k 



where note that = Vu. In other words, expanding (VQ + ^iVi +■ ■ ■ + V^))'^ generates 



all the possible vector fields Vw with w 
coefficients of powers of one-half. 



and vl(w) = k, with the appropriate 



Proof In the series solution for the flow-map ift, all stochastic information is encoded in 
the words on the left and the geometric information on the right. Taking the expectation 
of the flow-map, noting that expectation is a linear operator, and using the two results 
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above for E Jw and the combinatorial expansion, we get 
Eipt = E ^ JwVw 

loGA* 

k>On{w)=k 

k>0 
— exp(t£). 

Now note that exp(f/l) generates the semi-group for the solution to the parabolic 
differential equation in the theorem. □ 

Example (Heston model). In the Ifeston model [18], a stock price St is modelled by a 
stochastic process xt — log St with variance process vt which evolve according to: 

dxt = /idt + ^/iHdWl, 

dvt = K,{e - vt) dt + e ^t [pdWt^ + ^l-p^dW?) , 

given in Ito form. Note that the variance is a mean-reverting process; it tries to revert 
to the mean value 9 at rate k. Using the Feynman-Kac formula the corresponding 
partial differential equation for u{x,v,t) = E (^f{xt,vt) \ xq = x, vq = is 

Ut — p,Ux + n{6 — v)uv + ^vuxx + pevuxv + ^e'vuw 



Remark. The Feynman-Kac formula shows how we can solve a partial differential equa- 
tion to obtain information about the solution yt to the stochastic differential equation 
at time t £ [0, T]. In the reverse direction, to numerically solve high dimensional diffu- 
sion problems in the form of deterministic partial differential equations, we need only 
simulate an A''-dimensional stochastic equation for yt and then compute the expectation 
Eyt to find the solution. 



5 Monte-Carlo simulation 

In Monte-Carlo simulation, we generate a set of suitable multidimensional sample 
paths say W{ijj) ■— (W^{ijj), . . . ,W'^{ijj)) on [0,r]. We generate a large finite set of 
paths, each labelled by ui. Here the number of paths, say P, must be large enough 
so that, for example, any statistical information for the solution yt that we want to 
extract is sufficiently robust. For each sample path W{u)), we generate a sample path 
solution to the stochastic differential equation on [0, T]. This can often only be 
achieved approximately, by using a truncation of the stochastic Taylor expansion for 
the solution y on successive small subintervals of [0,r]. Having generated a set of 
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approximate solutions yt{i^i) at time t G [0, T] for every cji for i = 1, . . . , P, we estimate 
the expectation E/(j/t) by computing 

P 

i=l 

regarded as a suitable approximation for 

E/(yt) := J f{yt{oj))dP{uj) 

over all possible paths. Now a natural question arises. Do the suitable paths W{ljj) we 
generate, have to be sample Brownian paths to compute the mean above, or can we 
choose different paths that will still generate the expectation effectively? We discuss the 
latter case (weak simulation) briefly next. We then move onto a more comprehensive 
treatment of the former (strong simulation) case, which also allows us to include path- 
dependent features. 

5.1 Weak simulation 

There are several weak simulation strategies to approximate Ef{yt). The simplest 
and most common is to replace the driving Wiener paths by paths generated as 
follows. Construct paths by generating increments Z\VF*(tn, in+l) over the computation 
subinterval [tn,t„+i] by the binomial branching process 

where h = tn+i — tn- Then, depending on the ordinary numerical integration scheme 
employed for each such path, one can show that for some order of convergence p, for 
some t G [0, T], we have 

\\Efiyt)~Ef{yt)\\=OihP). 

See Kloeden and Platen [22] for more details. Another promising method is the cubature 
method of Lyons and Victoir [2H] • 

5.2 Strong simulation 

In a strong simulation, discrete increments AW^ in each computation interval are 
directly sampled from the Gaussian distribution. Indeed we generate each multidimen- 
sional path W{uj) by choosing increments in each component as follows 

AW\tn,tn+l) ~ ■ N(0,1). 

This is more expensive, but sample paths W{ijj) generated in this way, allow us to 
compare yt(a;) and i/t(aj) directly in the sense that one can show 

E\\yt-yt\\=0{h^) 

for some order of strong convergence p/2; which we will discuss in detail in Section [T] 
Often in practice we take || ■ || to be the Euclidean norm so that the convergence shown 
is in the L^-norm. 
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Given a sample mulitdimensional path W{uj) on [0, T], how do we actually construct 
an approximate solution yt? Here we are guided by the stochastic Taylor expansion. In- 
deed, classical strong numerical methods are based on truncating the stochastic Taylor 
expansion 

yt = ^ Jw{t) Vwiyo), 

uiGA* 

and applying the approximation over successive subintervals of the global interval of 
integration [0, T]; see Kloeden and Platen [22] or Milstein [31]. We present three simple 
example numerical approximation methods. 



5.3 Euler-Maruyama method 

If we truncate the Ito form of the stochastic Taylor series after the first order terms we 
generate the Euler-Maruyama numencal method as follows: 

d 

Vn+l = j/n + hVo{yn) + ^ (ZiVF*(tn, t„+l)) Vi{yn). 
i=l 

2_ 

This is a numerical scheme with global order of convergence of /i 2 . We explain in 
Section [3 why we have used the Ito drift vector field here. 



5.4 Milstein method 

We now truncate the stochastic Taylor series after the second order terms. This gen- 
erates the Milstem numencal method given by 

d d 

yn+l=yn + hVo{yn)+^^{AW\tn,tn+l))Vi{yn)+ ^ J^j (t„ , t„+i ) (j/„) . 
i—1 i,J — 1 

An important and expensive ingredient in this method, is the simulation of the multiple 
integrals Jij{tn,tn+i) for i ^ j shown, on each integration step. When i = j, the 
multiple integrals Jii{tn,tn+i) are cheaply evaluated by 

Mtn,tn+l) = J^"*' ' dW^.dW^, = ^{AW\tn,tn+l)'f. 

When j 7^ j we have by integration by parts that 

Hence we need only compute one double integral for each pair i ^ j. Equivalently we 
need only compute the Levy area given by 

since Jij = \JiJj + A^j. By Stokes' Theorem, the Levy area on the interval [tn,tn+i] 
is the chordal area for the path {W^, W-') on [tn,tn^i]. This can be observed directly 
by definition since 



Jo 



We consider the issue of simulating the Levy area in some detail in Section [B] The 
Milstein scheme has global order of convergence h. 
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Binomial path samples Brownian patin samples 




0.5 1 0.5 1 

time time 

Weal< solution samples and expectation Strong solution samples and expectation 




time time 



Fig. 2 Weak and strong simulations for the scalar linear example given in Section 2, with a = 3 
and 6 = 2. For this example we took yo = I, T = 1, h = 0.05 and P = 10 for both simulations, 
though only 5 sample paths are shown above in all cases. Top left are 5 sample binomial 
branching paths w, and top right are 5 sample Brownian paths W. Lower left are 5 sample 
solution paths y using the binomial branching paths, while lower right are 5 sample solution 
paths Y using the Brownian paths; both computed using the Euler-Maruyama method. At 
each time-step the thick black line shows the average value over P = 10 samples and the red 
line is the analytic solution for the expectation. 

5.5 Castell-Gaines method 

Consider the exponential Lie series ipt = logiySt generated by taking the logarithm of 
the stochastic Taylor series for the flow-map, i.e. 

tPt = (V5t - id) - ^{ft - idf + ^(<pt - idf + ■■■ 
d 

i—O i>j 

This series is also known as the Chen-Strichartz, Chen-Fleiss or Magnus series. The 
Castell-Gaines method is a strong numerical method based on truncating the expo- 
nential Lie series. As for the methods above, we generate a set of multidimensional 
paths W{ijj) on [0,r] with Wiener increments AW^{tn,tn+i) sampled on the scale 
h = — tn- On each computation interval we relace the Ji{tn,tn+i) by 

the Normal samples AW^{tn,tn+i)- If required, the Levy area increments Aij{tn,tn+i) 
shown, are also replaced by suitable samples Aij{tn,tn+i) as we outline in Section (6] 
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Then across the computation interval [tn,tn+i], we have 
d 

i—0 i>j 

The solution at time tn+i is then approximately given by 

J/t„+i « exp(V)t„,t„^i) oyt„. 

Note that for each path AW^{tn,tn+i) and Aij{tn,tn+i) axe fixed constants. Hence 
the truncated Lie series tpt„,t„+i is itself an autonomous vector field. Thus, for r G [0, 1] 
and with u{0) = yt^ , we solve the ordinary differential equation 

u'{r) = Tpt^j^+i o u{r)- 

Using a suitable high order ordinary differential integrator generates u{l) ~ ytn+i- 

Without the Levy area the Castell-Gaines method has global order of conver- 
gence /i2 ; while with the Levy area it has global order of convergence h. Castell and 
Gaines |10llllj prove that their strong order h'^ method is always more accurate than 
the Euler-Maruyama method. Indeed they prove that this method is asymptotically 
efficient in the sense of Newton [34]. Further in the case of a single driving Wiener 
process (d = 1), they prove the same is true for their strong order h method. By 
asymptotically efficient we mean, quoting from Newton, that they "minimize the lead- 
ing coefficient in the expansion of mean-square errors as power series in the sample 
step size". 

6 Simulating the Levy area 

A fundamental and crucial aspect to the implementation of strong order one or higher 
integrators for stochastic differential equations, is the need to successfully simulate the 
Levy chordal areas Aij{tn,tn+i), when the diffusion vector fields do not commute. 
This aspect is more than just a new additional concern once we step off the cliff edge 
of simple path increment approximations with frozen vector fields characterized by the 
Euler-Maruyama approximation. It also represents a substantial technical difficulty. 
Here we will outline several methods employed to simulate it sufficiently accurately; the 
important distinguishing criterion for the success of the method will be its asymptotic 
rate of convergence as /i — >■ 0. A sample survey of these methods can be found in 
Ryden and Wiktorsson |36) . Here we will focus on the case of two independent Wiener 
processes and and the requirement to simulate Ai2ih) := Ai2(tn,tn+i), given 
the Normal increments AW^{h) := AW'^ (tn,tn+i) and AW'^{h) ■= AW\tn,tn+i) 
across [tn,tn+l]- 

6.1 Simulating Normal random variables 

We will start with the question: what is the most efficient method for generating 
AW^{h) and AW'^{h)7 The simple and direct answer is to use the Matlab command 



sqrt (h) *randn 
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This command invokes an algorithm that has been scrupulously refined and adapted 
over the years. One of the simplest efficient earlier incarnations of this algorithm is 
Marsaglia's polar method [30) . The Box-Miiller method is also very simple but not 
quite as efficient; see Kloeden and Platen [22] for a discussion of these issues. Also 
see Moro's inversion method ^S]- We outline Marsaglia's method here because of its 
simplicity and effectiveness. 

Algorithm 1 (Marsaglia's method) To produce two standard Normal samples: 

1. Generate two independent uniform random samples Ui,U2 G Unif([— 1, 1]); 

2. If S ■= Ui + U2 < i continue, otherwise repeat Step 1; 

3. Compute Xj — Uij \J—2 ln(5')/S, for i — 1,2; then X\ and X2 are independent 
standard Normal samples. 

6.2 Conditional distribution of Levy area 

The characteristic function <j) of the probability density function for ^ = A12 (h) given 
AW^ih) and AW'^ih) is 

■^(^^ = -T7r^^^p(-^»'(i'^Ccoth(iftO- 1)) 

smh(^/!,4) V / 

where = [[AW'^{h)) + {AW^{h)) ) /h. Levy derived this in a very succinct calcu- 
lation in 1951; see Levy [241 pp. 171-3]. Since is the characteristic function, i.e. the 
Fourier transform of the corresponding probability density function, the actual prob- 
ability density function (f) is given by the inverse Fourier transform (see for example 
Gaines and Lyons [IS]*): 

/•OO 

= W <^(Ocos(a;e)dC. 
Jo 

The ungainly form of this probability density function means that generating samples 
is not likely to be easy. For example, the simplest method for sampling from a contin- 
uous distribution / is based on the inversion of its (cumulative) distribution function 
F{x) := J^QQ /(?)) dr]. If we sample from the uniform distribution, say U ~ Unif([0, 1]), 
then F~^{U) is a sample from the target distribution. For this to be a practical sam- 
pling method we must have an analytic form for F or an extremely efficient quadrature 
approximation for the integral in F at our disposal. We don't have this for the proba- 
bility density function of the Levy area (p. 

Several methods have been proposed for sampling from (j). Gaines and Lyons [13] 
proposed one of the most efficient, based on Marsaglia's rectangle-wedge-tail method. 
However it can be complicated to implement. Kloeden and Platen [22] and Wiktorsson 
have proposed methods based on the Karhunen-Loeve expansion are much easier to 
code. Ryden and Wiktorsson [30] proposed a method based on recognising the charac- 
teristic function as a product of characteristic functions for a logistic random variable 
and an infinite sum of Poisson mixed Laplace random variables. Gaines and Lyons [14] 
also proposed a method based on the conditional expectation of the Levy area, condi- 
tioned on intermediate Wiener increments. We discuss these methods in the following 
four sections. Stump and Hill 39J have also proposed a very efficient method, whose 
potential in a practical implementation is yet to be explored. 
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Fig. 3 Sample two-dimensional Wiener path and enclosed chordal Levy area. 
6.3 Karhunen-Loeve expansion method 

Levy derived the form for the characteristic function (j) for the Levy area, using the 
Karhunen-Loeve expansion for a Brownian bridge. This is an expansion in orthogonal 
basis functions. The details can be found in Levy [24] or Kloeden and Platen [22| . 
If (7^; , V^j , Xj, , Yj. are independent N(0, 1) samples, also independent of AW^{h) and 
AW'^{h), then the Levy area can be represented by 

oo 

fe=i 

In practice, we truncate this expansion to only include k ^ Q terms and use the 
truncation, Ai2{h), as an approximation for the Levy area. The important question 
now, is as far as strong simulation is concerned, how many standard Normal random 
variables do we need to simulate in order to have a sufficiently accurate Levy area 
sample Ai2{h), i.e. how large must Q be? Note that the coefficients in the above 
expansion scale like h/k. The properties of the tail of the series, i.e. all the terms for 
k ^ Q + 1, mean that it scales as hj^fQ. For a Milstein numerical approximation we 
require that the strong error is locally of order ft 2 ; see Section [T] Hence we must choose 
Q « for a sufficiently accurate sample. 



6.4 Ryden and Wiktorsson's method 

Ryden and Wiktorsson [36) proposed several methods, we detail here the most expe- 
dient. The characteristic function i\> is the product of two characteristic functions 

<^X(/i)(C) = — and = exp(^-ia^(i/iCcoth(i/iC) - 1)) 

corresponding to the random variables XiK) and Yijn), respectively. We observe that 
4>x(h) is the characteristic function of a logistic random variable which can be generated 
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by the inverse method, i.e. pick U ~ Unif([0, 1]) andlet X(h) = (/i/Stt) log ((7/(1 - (7)) . 
Then using the identity 



z coth z — 1 — 2 



z^ 



fc = l 



we observe that 

e 



^^^^ {27vk/hy+e 



This can be viewed as a sum of compound Poisson random variables. Indeed if for each 
fc G N, we generate A^j. ~ Poisson(a^) and then, for j = 1, . . . , A^j. generate independent 
Laplace random variables Yjf^ ~ Laplace(l/A;), then 

oo Nk 

k=ij=i 

has density (j>Y{h)- In ^ practical implementation we truncate this expansion to include 
k ^ Q terms, and use the truncation as an approximation for the Levy area. Further 
the tail sum, by the central limit theorem, is asymptotically Normally distributed and 
can be approximated by a Normal random variable. This provides quite a dramatic 
improvement as it is possible to show that this method only requires the number of 
standard Normal samples to be Q . 



6.5 Wiktorsson's method 

Wiktorsson proposed a method that uses the Karhunen-Loeve expansion method, but 
also simulates the tail sum as in the last method. Again, by the central limit theorem, 
the tail sum can be approximated by a Normal random variable, and the corresponding 
improvement is that this method only requires the number of standard Normal samples 
to be Q « Wiktorsson's method has been successfully implemented by Gilsing 

and Shardlow [16) in their SDELab, to where the interested reader is referred. 



6.6 Conditional expectation 

One more approach to simulating the Levy area, or equivalently Ji2{tn, tn+i), is based 
on replacing Ji2{tn,tn+i), by its conditional expectation Ji2{tn,tn+i), as follows. 
Suppose we are about to perform the numerical update for the solution across the 
interval [tn,tn+i]- We generate Q pairs of independent standard Normal random vari- 
ables Xq,Yq ~ N(0, 1) for g = 1, . . . , Q. Set = t„ + qAt, for q = 0, . . . , Q - 1, 
and AW^{Tq) = ^^AiXq and AW^{rq) = VAtVq, where At is defined by QAt = h. 
We thus generate a two-dimensional Brownian sample path on [tn,tn+i]- We can take 
AW^{h) and AW'^(h) to be the increments across the interval [tn,tn+i]- More impor- 
tantly, we can use the intervening path information we have generated on the scale At 
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to approximate Ji2{tn,tn+i)- Indeed, Ji2{tn,tn+i) can be expressed as 

Jl2itn,tn+l) = ^ ^ dW^, dW? 

= J2 - ) + Wl ) dW? 

q=0 -'^l 

0-1 Q-1 

The quantity 

0-1 

9=0 

represents the expectation of Ji2{tn,tn+i) conditioned on the increments AW^{Tq) 
and AW'^{Tq). Prom an algebraic and geometric perspective, Jn{tn-,tn+i) represents 
a suitable approximation to Ji2(tn, tn+i)- Computing its mean-square strong error we 

see that 

2 

\\Jl2{tn,tn+l) - Jl2{tn,tn+l)\\j^^ = ^ II '^2 (rg, T^+i) ||^, = = /l^Q- 

q={) 

Hence its root-mean-square strong error is hj^fQ. Thus, as for the Karhunen-Loeve 
expansion approach, to achieve a suitable approximate sample for the stochastic area 
integral, this method requires Q « hT^ . One advantage of this method is that it is 
very convenient for generating log-log error plots. 



7 Strong error 

We will focus here on the global, strong l? error. In practical terms, the global error 
is generated by the accumulation of contributions from the local error. The local error 
is itself the leading order terms in the remainder Rt, say, of our truncated stochastic 
Taylor scries. Note that there is also a contribution to the global error from the approx- 
imate simulation of the Levy area, however we will assume here, that the Levy area 
has been sufficiently accurately simulated, as discussed in detail in the last section, so 
that its contribution is small, in comparison to the truncation error. 

Suppose we base a strong numerical approximation on truncating the stochastic 
Taylor expansion (in Stratonovich form) . Let y denoted the truncated expansion and R 
the corresponding remainder; hence the exact solution is y = y + R. To guarantee our 
numerical scheme based on such a truncation is globally of order /i"*, where m € Z/2, 
which terms must we keep in y? We give the following rule. 

Rule of Thumb: Terms in the remainder R of L'^ measure /i™ with: 

— zero expectation, accumulate so as to contribute to the global error as h™~'^ order 
terms; 

— non-zero expectation, accumulate so as to contribute to the global error as h"^~^ 
order terms. 
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Hence to achieve an integrator with global error of order /i™, we must retain in y: 

— all terms with LF' measure of order /i™ for all m' ^ m; 

— the expectation of all terms of order /i^^s which have non-zero expectation (the 
corresponding terms left in remainder will then have zero expectation). 

Example (Euler-Maruyama). Recall that we based the Euler-Maruyama approximation 
on the truncated Ito Taylor series. If we had truncated the stochastic Taylor series in 
Stratonovich form, then according to the rules above, to achieve global order we 
should retain in our integrator the expectation of the terms 

d 

Y,{V^■^yV{)Ju■ 

i=l 

Since E Jjj = ^h, we thus recover the corresponding truncated Ito Taylor series. 
8 Further issues 

There are many important simulation issues we have not had space to discuss. Chief 
among these is the numerical stability of the strong methods we have explicitly outlined. 
This issue is discussed in Higham [TS] and more can be found for example in Buckwar, 
Horvath-Bokor and Winkler [S]. 

A Stratonovich to ltd relations 

We give here some Stratonovich to Ito relations for convenience for the reader — more details 
can be found in Kloeden and Platen 1221 . For the words w shown, the Stratonovich integrals 
Juj can be expressed in terms of Ito integrals 1^ as follows: 

W = aia2 ■■ Jw = Iw + \lo '5ai=a2#0; 
W = 010203 : = + |(/0a3 <5ai=a27^0 + -^aiO Sa2=a3jto)j 

W = 01020304 : Ju, = 7„ -I- i/oO ^ai=a2#0<5a3=a47^0 

+ ^(lOa^ai Sai=a2j^0 + -faiOa4 ^a^^asjtO + -faia20 ^a^^a^j^^o)- 

Note that the expectation of any Ito integral Im is zero, i.e.: Elm = for any word € A* 
which has at least one non-zero letter. 

B Sample program for weak and strong Euler-Maruyama 

We provide the listing for the weak vs strong Eulor— Maruyama simulation shown in Figure [2] 
Listing 1 Weak vs strong simulation 



XX Weak and strong simulation example plot 

XX Euler-Maruyama approximations 
XX Example scalar linear SDE 



XX 



Parameter values 
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a=3.0; % drift coefficient 

b=1.4; % diffusion coefficient is b 

yO=1.0; % initial data 

P=10; Z total # of sample paths 

h=0.05; % stepsize 

T=1.0; Z global time interval 

N=T/h; % number of subintervals 

XX Binomial branching process increments 

dw=zeros (N , P) ; 
w=zeros (N+1 , P) ; 

binom=binornd (1 , 1/2 , [N , P] ) ; X Gives or 1 with prob 1/2 

dw=sqrt (h) * (1-2* binom) ; X Binomial increments dw 

w (2 : N+1 , : )= cumsum (dw , 1) ; X Binomial paths themselves 

XX Weak solution by Euler -Maruyama approximation 

y=zeros (N+1 , P) ; X y is weak solution 

for p=l:P % set loop for each path 

y(l,p)=yO; X initial data 

for n=l : N 

y(n+l ,p)=y(n,p)+a*y(n,p)*h+b*y(n,p)*dw(n,p) ; 

end 

end 

%% Approximate Wiener path increments 

dW=zeros (N , P) ; 
W=zeros (N+1 , P) ; 

dW= sqrt (h ) * r andn (N , P ) ; X Brownian increments dW 

W (2 : N+1 ,:)= cumsum (dW , 1) ; X Brownian paths themselves 

XX Strong solution by Euler -Maruyama approximation 

Y=zeros (N+1 , P) ; X Y is strong solution 

for p=l:P X set loop for each path 

Y(l,p)=yO; X initial data 

for n=l : N 

Y(n+1 ,p)=Y(n,p)+a*Y(n,p)*h+b*Y(n,p)*dW(n,p) ; 

end 
end 

XX Compute the expectations of y and Y at each timestep 

expect_y=zeros(N+l,l); 
expect _Y = zeros (N + 1 , 1) ; 

for n=l:N+l 

expect_y (n)=mean(y(n, :)) ; 
expect_Y(n)=mean(Y(n, :)) ; 

end 
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Fig. 4 Error vs stepsize and error vs CPU time for the Heston model. The parameter values 
can be seen in the program listing. 



C Example strong simulation program 
C.l Heston model strong simulation 

We provide here a sample program that shows how to perform log-log error plots for a strong 
simulation. We used a real example, the Heston model, and applied the full truncation Euler— 
Maruyama type numerical scheme devised by Lord, Koekkoek and Van Dijk 1251 . The log-log 
error vs stepsize, and error vs CPU time, are shown in Fig. [l] Note that to estimate the strong 
global error, we must compare solutions for different stepsizes along the same path, before 
taking the expectation. 

Listing 2 Heston model strong simulation 



%% Heston model integrati 

'/, Global data 

T0 = 0; 
T = l ; 
M = 10; 

hmin=(T-T0)/2-M; 
Mstart =4; 

hmax=(T-T0)/2"Mstart ; 
Q=(l/(hmin))*(hmax/hmin) ; 
dt=hmax/Q ; 
R = M-Mstart +1 ; 

P=100 ; 

alpha=2 . ; 

theta=0 . 09 ; 

beta=0 . 1 ; 

rho =0.5; 

mu=0 . 05 ; 

ic = [1 . ; . 09] ; 

YFT = zeros (P ,R ,2) ; 
clockYFT=zeros (1 , R) ; 



% integration interval... 
% . . . is [TO , T] 

% smallest timestep 

% largest timestep 

Z ^ of quad pts for hmin 

X quad scale for hmin 

% # of solution approximations .. . 

% . . . at different stepsizes 

X total # of Brownian paths 

X parameters 



% initial data 

X approximate solution 
X CPU timings 
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for p=l:P 

for r=l:R 

YFT (p , r , : ) = ic ; % make sure start with IC 

end 

end 

for j j =1 : 2" Mstart 
YFTold=YFT ; 

for p=l:P % loop for each path 

X'/o Start loop computing over intervals of length hmax 

siv=hmax/hmin ; 

%% Generate dWl and dW2 on smallest scale 

dWl = sqrt (dt ) *randn (1 , Q) ; '/, Brownian increments dWl 

dW2 = sqrt (dt ) *randn (1 , Q) ; '/, Brownian increments dW2 

dWO=( zeros (1 , Q)+l)*dt ; 

'/,'/, Start loop computing Js with different stepsizes 



for r=l:R 
SF=2- (r-1) ; 
h=SF*hinin ; 
L=hmax/h ; 
QR=Q/(SF-2) ; 

% Compute dwO , dwl , dw2 , wO , ml, 

dwO=zeros (1 , QR) ; 
dwl=zeros (1 , QR) ; 
dw2=zeros (1 , QR) ; 
wO=zeros (1 , QR) ; 
wl=zeros (1 , QR) ; 
w2=zeros (1 , QR) ; 



'/, scale factor for stepsizes... 
'/, stepsize 

'/, # of timesteps needed 

% ^ of quadrature steps (total)... 

w2 



for j=l:QR 

dwO ( j ) = sum (dWO ((j-l)*(SF'"2) + l:j*(SF-2))); 
dwl (j)=sum(dWl ((j-l)*(SF-2)+l:j*(SF-2))); 
dw2 (j)=sum(dW2 ((j -l)*(SF-2)+l: j*(SF-2))) ; 

end 

QF = qR/L; '/, # of quadrature steps in h... 

for n=l : L 

wO((n-l)*QF + l:n*qF) = cumsum ( dwO ((n-l)*QF + l:n*QF)) ; 
wl((n-l)*QF+l: n*qF)= cumsum (dwl((n-l)*QF+l:n*QF)); 
w2((n-l)*QF+l: n*QF)= cumsum (dw2 ((n-l)*QF+l:n*qF)); 

end 

oldclockYFT = clockYFT (r) ; 
ts=cputime ; 

YFT (p , r , : ) = YFTappr ox (p , h , L , TO , QF , dwO , dwl , dw2 , wO , wl , w2 , ... 

alpha , theta , beta , rho ,mu , . . . 

YFTold(p,r , :)) ; 
clockYFT (r) = cputime -ts + oldclockYFT ; 



end 
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end 
end 

stepsizes=loglO((2."([l:R-l])) *hmin) ; 

saveC'stepsizes ' , 'stepsizes ') 

save ( ' P ' , ' P ' ) 

save ( ' YFT ' , ' YFT ' ) 

save ( ' clockYFT ' , ' clockYFT ' ) 



G.2 Program listing: integrator 

We give the program listing for the Heston model full truncation Euler-Maruyama integrator. 
Listing 3 Pull truncation Euler-Maruyama integrator 



Z% Full truncation for volatility and exponential for index: 

function trunc = YFTapprox (p , h , L , TO , QR , dWO , dWl , dW2 , . . . 

WO , Wl , W2 , alpha , theta , beta , rho ,mu , ic ) 

t rune = zeros (2,1); 

% Compute dJl , dJ2 on scale h=SF*hmin 

Jl=zeros (1 , L) ; 
J2=zeros (1 , L) ; 
Jl (1)=W1 (QR) ; 
J2 (1)=W2 (QR) ; 
pts=2*C)R: QR:L*qR; 
Jl (2: L)=W1 (pts) ; 
J2(2:L)=W2(pts); 

%% initial data 

S=ic(l); % asset price 

v=ic (2) ; % volatility 

for n=l : L 
Sold=S; 
vold=v ; 

S=exp ( (mu-max (0 , void) /2) *h+sqrt (max (0 , void ) ) * Jl (n) ) * Sold ; 
v=vold+ alpha* (theta -max (0,vold))*h+beta*( rho * Jl (n) . . . 
+ sqrt (1-rho "2) * J2 (n) )* sqrt (max (0 , void) ) ; 

end 

trunc = [S ; v] ; 
end 



C.3 Program listing: log-log strong error plots 

The following program performs the log-log plots for the strong error measure. 
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Listing 4 Log-log strong error plots 



%X Log-log strong error plots 

load stepsizes 
load YFT 
load clockYFT 
load P 

R=length( stepsizes ) ; 
er r or YFT = zeros (1 , R) ; 
diffYFT=zeros(P,R); 

% the "norm" below is the Euclidean norm 

for r=l:R 

for p=l:P 

dlffYFT(p,r)=norm(YFT(p,r+l)-YFT(p,l)) ; 
end 

end 

% now take the L'2 norm measure 
errorYFT=sqrt (mean (diff YFT ."2,1)); 
f igure 

subplot (1,2,1) 

plot ( stepsizes (1 : end-1) , loglO ( errorYFT (1 : end -1) ) , . . . 

'-ks' , 'LineWidth ' ,2) 
xlabel('log_{10}(stepsize) ') 
y label ( ' log. {10} ( global u err or ) ' ) 
t it le ( [ ' Number u of u sampled upaths = ',int2str(P)]) 
subplot (1,2,2) 

plot (loglO (clockYFT (2: end-1) ) , loglO ( errorYFT (1 : end-1) ) , . . . 

' -ks ' , ' LineWidth ' ,2) 
xlabel ( ' log_{10} ( CPUutime) ' ) 
ylabel (' log_ {10} ( global u error )' ) 
t it le ( [ ' Number u of u sampled upaths = ',int2str(P)]) 
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